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A minimal kinetic model is used to study analytically and numerically flows at a mi- 
crometer scale. Using the lid-driven microcavity as an illustrative example, the interplay 
between kinetics and hydrodynamics is quantitatively visualized. The validity of various 
theories of non-equilibrium thermodynamics of flowing systems is tested in this nontrivial 
microflow. 



1. Introduction 

Flo ws in microdevices constitute an emerging application field of fluid dynamics l|Beskok fc Karniadakisl 
l200l|) . Despite impressive experimental progress, theoretical understanding of microflows 
remains incomplete. For example, even though microflows are highly subs onic, the as- 
sumpt ion of incompressible fluid motion is often questionable (see, e. g. IZheng et all 
l)2002|) '). Interactions between the relaxation of density variations, the rarefaction and 
the flow geometry are not completely understood. One of the reasons for this is the lack 
of commonly accepted models for efficient simulations of microflows. 

In principle, microflows can be studied usin g^ molecular level methods, such as the Di- 
rect Simulation Monte Carlo (DSMC) method ()Birdll994|) . However, molecular dynamics 
methods face severe limitations for subsonic flows. The number of particles required for 
simulations in realistic geometries with large aspect ratios and th e number of time steps 
needed to reach the statistical steady state are prohibitively high ijOran et all\l99c^i (the 
time step of DSMC is ~ 10~ 10 sec, while relevant physical time-scales are ~ 10 — lOOsec). 
Moreover, a detailed microscopic description in terms of the particle distribution function 
is not necessary for design/endinccring purposes. Thus, development of reduced models 
enabling efficient simulations is an important issue. 

In this paper we consider a two-dimensional min imal kinetic model with nine discrete 
velociti es (IKarlin et q/Jl999t I Ansumali et a?j2003l). It has been recently shown by several 

ill 



group s l|Ansumali fc Karlinl2002 fell Ansumali et a?j2004tlNiu eraJ2003lSucci fc Sbragariial 
2004 ) that this model compares well with analytical results of kinetic theory ijCercignanl 



1975() in simple flow geometries (channel flows), as well as with molecular dynamics 
simulations. The focus of this paper is the validation of several theoretical models of 
(extended) hydrodynamics versus direct numerical simulation in a non-trivial flow. For 
that purpose, we considered the flow in a lid-driven cavity. We will visualize the onset of 
the hydrodynamic description, the effect of the boundaries etc. 

The paper is organized a s follows: For completeness, the kinetic model ()Karlin et all 
Il999t lAnsumali et all\2()()3) is briefly presented in sectional In sectional we show the 
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relation of our mo del to the we ll-known Grad moment system derived from the Boltzmann 
kinetic equation <jGradlll94flh . We compare analytically the dispersion relation for the 
present model and the Grad moment system. This comparison reveals that the kinetic 
model of section [21 is a superset of Grad's moment system. In section 01 a parametric 
numerical study of the flow in a micro-cavity is presented. Results are also compared to 
a DSMC simulation. In section the reduced description of the model kinetic equation 
is investigated, and a visual representation quantifying the onset and breakdown of the 
hydrodynamics is discussed. We conclude in section wit h a spectral analysis of the 
steady state flow and some suggest ions for further research ijTheodoropoulos et aZJl200ol 
l2004HKevrekidis et al.ll2003L 120041) . 



2. Minimal kinetic model 

We consider the discrete velocity model with the following set of nine discrete velocities: 



= [0,1, 0,-1, 0,1 -1,-1,1] 



= [0,0,1,0,-1,1,1,-1,-1]. 



(2.1) 



The local hydrodynamic fields are defined in terms of the discrete population, f iy as: 



^ ] c xi: Cyi\ — {p, j Xl Jj,}, 



(2.2) 



i=i 



where p is the local mass density, and j a is the local momentum density of the model. 
The populations /, = f(x,Ci,t) are functions of the discrete velocity Cj, position x and 
time t. We consider the following kinetic equation for the populations (the Bhatnagar- 
Gross-Krook single relaxation time model): 



dtfi + Ci ■ d x fi 



1 



(fi-fP(f)), 



(2.3) 



where r is the relaxation time, and /? q is the local equilibrium l|Ansumali et al.ll2003|) : 



f p - m ( 2 - yrw) ( 2 - ^) (^t5Z) (^^Z 

(2.4) 

where u a = j a /p, and the speed of sound is c s = 1 / V3. The local equilibrium distributio n 
/j 6q is the minimizer of the discrete H function <|Karlin 



A 



with weights W = 



16 4 4 4 4 1 1 1 1 



36' 36' 36' 36' 36' 36' 36' 36' 36 



, (2.5) 



under the constraints of the local hydrodynamic fields (|2.2(l . Note the important factor- 
ization over spatial components of the equilibrium l|2.4(l . This is similar to the famil- 
iar property of the local Maxwell distribution, and it distinguishes (|2.4J) among other 
discrete-velocity equilibria. 

In the hydrodynamic regime, the model recovers the Navier-Stokes equation with vis- 
cosity coefficien t u = yr, where y = pet is the pressure (see section EJ. The diffusive- wall 



approximation I Ansumali & Karlinll20025l) is used for wall boundary conditions 
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3. Grad's moment system and the minimal kinetic model 

3.1. The moment system 

It is useful to represent the discrete velocity model (|2.3[1 in the form of a moment system. 
For simplicity, we shall consider the linearized version of the model. Note that lineariza- 
tion is required only for the collision term o n the right h and side of equation H2.3|) . This 
is at variance with Grad's moment systems ijGradlll Q4PJ) where the advcction terms are 
nonlinear as well. We choose the following nine non-dimensional moments as independent 
variables: 



M 



P_ 

Po' 



Jx 
PoCs 



where 



Jy_ 

' PoCs 

ib = R 



P 



N P x 



<l.r 



4' 



Po<% ' pool p cl ' 2p cl ' 2p cf ' 2p cf 



(3.1) 



yyyy 



2R 



(3.2) 

is a scalar obtained from 4 t/l -order moments Rap^e = $3*=i fi c aiCpiGyiC$i, N = Y^=i fi( c xi 
P yy )/2 is the difference of the normal stresses, P = X^ =1 /jcf is the trace 



xxyy, 
9 



i)/2 = {Px 



of the pressure tensor, and q a = X)j=i fi c aicf is the energy flux obtained by contrac- 
tion of the third-order moment Q a /3-y = Y^=i fi c ai c /3i c -yi- Time and space are made 
non-dimensional in such a way that for a fixed system size L they are measured in the 
units of mean free time and mean free path, respectively: x' = x/(LKn),t' = t/r, where 
Kn = TCs/L is the Knudsen number. The linearized equations for the moments M 13.1fl 
read (from now on we use the same notation for the non-dimensional variables) : 



d t P + dxjx + Oyjy = 0, 

d t jx + d x (P + N) + d y P xy = 0, 

fit jy + dxPxy + d y (P~N) = 0, 

d t P + d x q x + d y q y = (p- P), 
d t N + d x (q x - Q xyy ) - d y (q y - Q yxx ) = —N, (3.3) 

&tPxy ~\~ &xQyxx H~ QyQyyx Pxyi 
&tQx ~\~ ^x^xxaa H~ $yRxyaa — i^jx Qx) j 
&tQy ~\~ dx^xyaa H~ dyRyyaa " (P'Jy Qy) : 

dtip + d x (j x - q x ) + dy (j y - q y ) = (2p - ip) ■ 

Furthermore, by construction of the discrete velocities (|2.1(l . the following relations are 
satisfied: 

Qxyy — 2q x 3j x , Qy XX — ^Qy ^3yi (^•^ : ) 

R X yaa = 3P XJ/ , R xxaa = 3 (p + ±n\ - \b, R yyaa = 3 (p - ^Nj - \b. (3.5) 

Apart from the lack of conservation of the energy and linearity of the advection, equa- 
tion (|3.3|) is quite similar to Grad's two-dimensional 8-momcnt system f . However, in the 
present case a particular component of the 4 t/l -order moment is also included as a vari- 
able. In other words, Grad's non-linear closure for the 4* ,l -order moment is replaced by 

f The variables used in the D- dimensional Grad's system are density, D components of the 
momentum flux, D(D+ 1)/2 components of the pressure tensor and D components of the energy 
flux. The number of fields in Grad's system is 8 for D = 2 and 13 for D = 3. 
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an evolution equation with a linear advection term. We note here that while the formu- 
lation of boundary conditions for Grad's moment system remains an open problem, the 
boundary conditions for the extended moment system (system 13.31 13.41 and 13. 511 arc well 
established through its discrete- velocity representation i|2.3[) l|Ansumali fc Karlinl20026j) . 
We also note that like any other Grad's syste m the present mode l reproduces the Navier - 
Stokes equation in the hydrodynamic limit llKarlin et aZ.lfl99 1 1 Ansumali eTal\\200^\ . 
The moment system <|3.3[1 reveals the meaning of the densities appearing in model: The 
dimensionless density is the dimensionless pressure of the real fluid in the low Mach 
number limit, while the momentum flux density should be identified with the velocity in 
the incompressible limit. With this identification, we shall compare the moment system 
(|3.3|) with Grad's system. 



3.2. One- dimensional Grad's moment system 

Since energy is not conserved by the model (|2.3|) . the comparison will be with another 
Grad moment system which (for D = 3) is usually referred to as the 10-moment system 
t. For one-dimensional flows, the linearized Gra d's 10-moment system can be written as 
l|Gorban fc Karlinll2005l iKarlin fc Gorbar]l2002t : 

d t p + jd x u x = 0, d t u x + d x P xx = 0, d t P xx + 3d x u x = - (P xx - p) , (3.6) 

where 7 is the ratio of the specific heats of the fluid, and 7 = (D + 2)/D for a D- 
dimcnsional dilute gas. This model can be described in terms of its dispersion relation, 
which upon substitution of the solution in the form ~ exp (u>t + ikx) reads: 

cj 3 + uj 2 + 3fc 2 u + jk 2 = 0. (3.7) 

The low wave-number asymptotic represents the large-scale dynamics (hydrodynamic 
scale of Kn -c 1), while the high- wave number limit represent the molecular scales 
quantified by Kn S> 1. The low wave number (Kn <ti 1) asymptotic, oj\, and the large 
wave number (Kn ^> 1) asymptotic, a>h, are: 

= {^^fc 2 ± V7fc, -1 - (-3 + l)k 2 ] , c h = {tl+jl ± lV 3k, -2} . 

(3.8) 

The two complex conjugate modes (acoustic modes) of the 0(k 2 ) dynamics, are given by 
the first two roots of u>\, and represent the hydrodynamic limit (the Navier-Stokes approx- 
imation) of the model. The third root in this limit is real and negative, corresponding 
to the relaxational behavior of the non- hydrodynamic variable (stress): the dominant 
contribution (equal to —1) is the relaxation rate towards the equilibrium value, while the 
next-order correction suggests slaving of viscous forces, which amounts to the constitu- 
tive relation for stress ((— 3 + 7)/2fc 2 ). Furthermore, the k 2 dependence of the relaxation 
term justifies the assumption of scale separation (the higher the wave-number, the faster 
the relaxation). The real part of the high wave-number solution Wh is independent of 
k, which shows that the relaxation at very hi gh Knudsen number is the same for a ll 
wavenumbers (so-called "Rosenau saturation" l|Gorban fc Karlir]ll996USlemrodlll998h ). 
Thus, the assumption of scale-separation is not valid for high Knudsen number dynamics. 

f The variables used in this D-dimensional Grad's system are density, D components of the 
momentum flux, and D{D + l)/2 components of the pressure tensor, resulting in 6 and 10 
variables for D — 2 and D = 3, respectively. 



ity 5 

Op^ i i i i , i , i i i , 




to,. Rt[t», ,] 



1 2 3 4 5 6 

k 

Figure 1. Real part of the solutions of the dispersion relation (equation Roots u>2,3 and 

coi correspond to Grad's subsystem (equation 13. Gt ). The real- valued root ujq and the complex 
conjugate roots 0^2,3 are extended hydrodynamic modes. 

3.3. Dispersion relation for the moment system 

The dispersion relation for the one-dimensional version of the moment system H3.3(l (i.e. 
neglecting all derivatives in the y-direction) reads: 

(cj 3 + u? + 3k 2 uj + k 2 )(u> 3 + 2uo 2 + (3k 2 + l)u> + k 2 )(l + w)((l + uj 2 ) + 2k 2 ) = 0. (3.9) 

The real parts of the roots of this equation (attenuation rates Rc[w(fc)]) are plotted in 
Fig- Has functions of the wave vector k. It is clear that for one-dimensional flows, the 
dynamics of three of the moments (p, j x , and P) are decoupled from the rest of the 
variables, and follows of the dynamics of the one-dimensional Grad's moment system 
|3~B|) with 7=1. 

The similarity between Grad's moment system and the present model is an important 
fingerprint of the kinetic nature of the latter f. Note that in the case of two-dimensional 
flows, the agreement between the present model and Grad's system is only qualitative. 
The present moment system is isotropic only up to 0(k 2 ). Thus, the dispersion relation 
of the model 1)2.3(1 is expected to match the one of Grad's system only up to the same 
order. In the hydrodynamic and slip-flow regime addressed below, this order of isotropy 
is sufficient. In the presence of boundaries and/or non-linearities, it is necessary to resort 
to numerics. Below we use the entropic lattice Boltzmann discretization method (ELBM) 
of the model (|2~3|) . 



4. Flow in a lid-driven micro-cavity 

The two-dimensional flow in a lid-driven cavity was simulated with ELBM over a range 
of Knudsen numbers defined as Kn = Ma/Re. In the simulations, the Mach number was 
fixed at Ma = 0.01 and the Reynolds number, Re, was varied. Initially, the fluid in the 
cavity is at rest and the upper wall of the domain is impulsively set to motion with 
mid = c s Ma. Diffusive boundary conditions arc imposed at the walls ( Ansumali fc Karlinl 
1200261) . and the domain was discrctized using 151 points in each spatial direction. Time 

f Grad |Grad 1949j) already mentioned that moment systems are particularly well suited for 
low Mach number flows. Qualitatively, this is explained as follows: when expansion in the Mach 
number around the no-flow state is addressed, the first nonlinear terms in the advection are of 
order u 2 /c 2 ~ Ma 2 . On the other hand, the same order in Ma terms in the relaxation contribute 
ti 2 /(rCs) ~ Ma 2 /Kn. Thus, if Knudsen number is also small, nonlinear terms in the advection 
can be neglected while the nonlinearity in the relaxation should be kept. That is why the model 
12.31 - linear in the advection and nonlinear in the relaxation - belongs to the same domain of 
validity as Grad's moment systems for subsonic flows. 
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Figu re 2. Flow in a micro-cavity for Kn = 0.1 and Ma = 0.14: DSMC simulation ij.liang et al\ 
2003) (left) , velocity vector plot and density isolines from ELBM (solid lines) with the DSMC 
density isolines (dashed lines) superimposed (right). 

integration is continued till the steady state is reached; matrix-free, Newton-Krylov fixed 
point a lgorithms for the accelerated computation of the steady state are also being ex- 
plored l(TheodoroDOulos et qiJl20f)0T) . 

4.1. Validation with DSMC simulation of the micro-cavity 

In the hydrodynami c regime, the model was va lidated using results available from con- 
tinuum simulations l|Ansumali fc Karlir]l2002a|). Far hig her Kn ~ 0.1, we compared our 
results with the DSMC simulation of l|jiang et aZjl2003|) . Good agreement between the 
DSMC simulation and the ELBM results can be seen in Fig. [3 It can be concluded, that 
even for finite Knudsen number, the present model provides semi-quantitative agreement, 
as far as the flow profile is concerned. We remind here again, that the dimensionless den- 
sity in the present model corresponds to the dimensionless pressure of a real fluid so that, 
for quantitative comparison, the density of ELBM model should be compared with the 
pressure computed from DSMC. 

4.2. Parametric study of the flow in the micro-cavity 

Fig. |3| shows the dimensionless density profiles with the streamlines superimposed for 
Kn = 0.001, 0.01, 0.1. For Kn = 0.001 (Rc = 10), the behavior expected from continuum 
simulations with a large central vortex and two smaller recirculation zones close to the 
lower corners can be observed. As the Knudsen number is increased, the lower corner 
vortices shrink and eventually disappear and the streamlines tend to align themselves 
with the walls. 

The density profiles, as a function of Kn, demonstrate that the assumption of incom- 
prcssibility is well justified only in the continuum regime, where the density is essentially 
constant away from the corners. This observation is consistent with the conjecture that 
incompressibility requires smallness of the Mach as well as of the Knudsen number. In 
hydrodynamic theory, the density waves decay exponentially fast (with the rate of re- 
laxation proportional to Kn) leading effectively to incompressibility. Thus, it is expected 
that the onset of incompressibility will be delayed as the Knudsen number increases. 

5. Reduced description of the flow 

The data from the direct simulation of the present kinetic model were used to validate 
the effectiveness of different closure approximations of kinetic theory in the presence of 
kinetic boundary layers in a fairly non-trivial flow, and to gain some insight about the 
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(a) (b) (c) 

Figure 3. Density isocontours for (a) Kn = 0.001, (b)Kn = 0.01 , and (c) Kn = 0.1 (the 
variation of the density is 0.995 < p < 1.005). Superimposed are the streamlines. 

required modification of the closure approximations in the presence of boundary layers. 
In this section, we will present such an analysis for two widely used closure methods, the 
Navier-Stokes approximation of the Chapman-Enskog expansion and Grad's moment 
closure. 

5.1. The Navier-Stokes approximation 

The Chapman-Enskog analysis JChaoman & Cowlinell97fl of the model kinetic equation 
leads to a closure relation for the non-equilibrium part of the pressure tensor as (the 
Navier-Stokes approximation) : 

o- xy = -TC 2 s {d y j x + d x j v ). (5.1) 

Fig. 01 shows a scatter plot of the xy component of the non-equilibrium part of the 
pressure tensor P xy — P X y, versus that computed from the Navier-Stokes approximation 
of the Chapman-Enskog expansion 15.1(1 . The upper row is the scatter plot for all points 
in the computational domain, while the lower row is the scatter plot obtained after 
removal of the boundary layers close to the four walls of the cavity, corresponding to 
approximately 10 mean-free paths. In all plots, the dashed straight line of slope equal to 
one corresponds to Navier-Stokes behavior. These plots clearly reveal that the Navier- 
Stokes description is valid away from the walls in the continuum as well as in the slip flow 
regime. On the other hand, it fails to represent hydrodynamics in the kinetic boundary 
layer. 

Perhaps the most interesting observation from Fig. 01 is the coherent, curve-like struc- 
ture of the off-Navicr-Stokes points. These trajectories tend to the Navier-Stokes ap- 
proximation as to an attractive sub-ma nifold. This structure resembles the traject ories 
of solutions to the invariance equation l|Gorban fc Karlir]l2005l iGorban et a^jEoO^l ob- 
served, in particul ar, in a similar problem of derivation of hydrodynamics from Grad's 
systems (see, e. g. iKarlin fc Corhaul fetinfl . p. 831, Fig. 12). A link betwee n solutions 
to invariance equations and the present simulations is an intriguing subject for further 
study. In the next subsection, we shall explore Grad's closure approximation for the 
present flow. 

5.2. Grad's approximation 

In contrast to the Chapman-Enskog method, the Grad method has an advantage that 
the approximations arc local in space, albeit with an increased number of fields. As the 
analysis of section|3|suggests, the dynamics of the density, momentum and pressure tensor 
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Figure 4. Scatter plot of the non-equilibrium part of the off-diagonal component of the pres- 
sure tensor P xy — P°JJ and corresponding value computed from Navier-Stokes approximation 
<J xy 15. II for all points in the domain ((a) and (b)), and after the removal of the boundary 
layer corresponding to approximately 10 mean- free path ((c) and (d)). Fig. (a,c) correspond to 
Kn = 0.001, while Fig. (b,d) correspond to Kn = 0.01. Navier-Stokes behavior is indicated by 
the straight line of slope equal to one. 

are almost decoupled from the rest of the moments, at least away from boundaries. This 
motivates the Grad-likc approximation for the populations, 



Grad 



+ TTJ ( P aP ~ SaPPcl) (c la C l0 - C 2 s 5 a[} ) 



(5.2) 



The set of populations parameterized by the values of the density, momentum and pres- 
sure tensor (|5.2|) is a sub- manifold in the phase space of the s ystem (12.31 . and can be 
derived in a standa rd way using quasi-equilibrium procedures ijGorban fc Karlml 120051: 
iGorban et aZjEoO^) . After taking into account the time discretization, we find the closure 
relation for the energy flux: 



_Grad 



P_ 

2(i 



Jo. 



2^ a 



(5.3) 



In Fig. [S] the scatter plot of the computed energy flux q x and the discrete Grad's closure 
15.311 is presented. Same as in Fig.QJ the off-closure points in Fig. [S]arc associated 



-,Grad 



with the boundary layers. The comparison of the quality with which the closure relations 
arc fulfilled in Fig. ^ and Fig. [S] clearly indicates the advantage of that a Grad's closure. 
Various strategies of a domain decomposition for a reduced simulation can be devised 
based on this observation. A general conclusion is that for slow flo ws Grad's closure in 
the bu lk along with the discretization of the boundary condition l|Ansumali fc Karlinl 
12002^ is the optimal strategy for the simulation of microflows. 



6. Conclusions and further studies 

We considered a specific example of a minimal kinetic model for studies of microflows, 
and compared it with other theories of noncquilibrium thermodynamics in a nontriv- 
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Figure 5. Scatter plot of the computed energy flux, q x , versus Grad's closure, q x ra : (a) 

Kn = 0.001, (b) Kn = 0.01. 



0.9S 0.99 




Figure 6. (a) Leading eigenvalues of the minimal kinetic model at steady state (square: 
Kn=10" 4 , circle: Kn=10~ 3 , X: Kn=10~ 2 , +: Kn^lO" 1 ); (b) Scatter plot as in Fig. E^a) for 
a state perturbed away from the steady state along the leading eigenvector (Kn=10~ 3 ). 



ial flow situation. The close relationship between Grad's moment systems and minimal 
kinetic models was highlighted. For the case of a driven cavity flow, different closure ap- 
proximations were tested against the direct simulation data, clearly showing the failure 
of the closures near the boundaries. Grad's closure for the minimal model was found to 
perform better than the Navier-Stokes approximation. This finding can be used to reduce 
the memory requirement in simulations, while preserving the advantage of locality. 

fn this paper, we have explored two classical closures of kinetic theory. In the future we 
are going to consider other closures such as the invariance correction to Grad's clo sure, 
and especially closures based on spectral decomposition iGorhan fc KarlirJ [20051) . To 
that end, the ELBM code was coupled with ARPACK i Lehouca et aZJ Il998 ) in order 
to compute the leading eigenvalues and the corresponding eigenvectors of the Jacobian 
field of the corresponding map at the steady state. In all cases, the eigenvalues are 
within the unit circle (Fig.[B{a)). The leading eigenvalue is always equal to one (reflecting 
mass conservation), and the corresponding eigenvector captures most of the structure of 
the steady state. As the Knudsen number decreases, eigenvalues tend to get clustered 
close to the unit circle. This happens because when the Knudsen number is small the 
incompressibility assumption is a good approximation, and mass is also conserved locally. 
The very close similarity between fig. Efb) and fig. Ufa), reveals that states perturbed 
away from the steady state along the leading eigenvector are also described well by the 
Navier-Stokes closure. 
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